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The slider-block Burridge-Knopoff model with the Coulomb friction law is studied as an excitable 
medium. It is shown that in the continuum limit the system admits solutions in the form of the 
self-sustained shock waves traveling with constant speed which depends only on the amount of the 
accumulated stress in front of the wave. For a wide class of initial conditions the behavior of the 
system is determined by these shock waves and the dynamics of the system can be expressed in 
terms of their motion. The solutions in the form of the periodic wave trains and sources of counter- 
propagating waves are analyzed. It is argued that depending on the initial conditions the system 
will either tend to synchronize or exhibit chaotic spatiotemporal behavior. 

PACS Number(s): 83.50. Tq, 91.30.Mv, 83.50.By, 03.40.Kf 
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I. INTRODUCTION 

Propagating self-sustained waves (autowaves) and 
more complex spatiotemporal patterns are characteristic 
of excitable media of different nature. A typical exam- 
ple of such a phenomenon is burning of black powder 
in a safety fuse. When the fuse is ignited at one end, 
the exothermic reaction releases the heat which is then 
spread out by heat diffusion. Thus, the neighboring re- 
gions of the fuse ignite leading to self-sustained propaga- 
tion of the combustion front. The phenomenon of non- 
attenuated propagation of waves is in fact common for 
a variety of physical, chemical, and biological systems 
. Traveling waves are experimentally observed in 
semiconductor and gas plasma, semiconductor and super- 
conductor structures, combustion systems, active optical 
media, magnetic media under illumination, autocatalytic 
chemical reactions, nerve and heart tissue (see and 
references therein). 

In order for self-sustained waves to be feasible, the sys- 
tem must possess two basic ingredients. First, the sys- 
tem must be excitable, that is, there has to be a thresh- 
old below which the perturbation of the steady homo- 
geneous state of the system decays, while perturbations 
of larger amplitude grow. In the example above, a suf- 
ficient amount of heat is needed to ignite black powder. 
Second, there has to be coupling between the regions 
of the system at different points in space. In the case 
of the safety fuse such a coupling is provided by heat 
diffusion leading to spread of the temperature and ig- 
nition of powder in front of the combustion zone. Thus, 
the prototype systems exhibiting self-sustained waves are 
reaction-diffusion systems [ . 

Recently, it was pointed out that an entirely differ- 
ent class of systems may be considered as excitable [Q. 
These are elastic media with friction exhibiting stick-slip 
motion. Both experimental observations and numerical 
simulations show that such systems are capable of sup- 



porting steadily propagating solitary waves in the form 
of shocks [0-|o|. These systems are also of special in- 
terest because they are used for modeling the dynamics 
of earthquakes ]7j-|l7|] . It is clear that systems exhibiting 
stick-slip motion have both necessary ingredients of ex- 
citable systems. The threshold behavior here is due to 
static friction, which prevents any motion in the system 
until some critical amount of stress is accumulated. The 
coupling of the elements of the system at different points 
in space is due to the non-locality of elastic stress. 

Singular perturbation techniques proved to be very ef- 
fective in treating problems of traveling wave propagation 
in reaction-diffusion systems p|-|^, |l8| , |l9| . These methods 
use strong separation of time scales in the problem to de- 
compose the dynamics of the system into fast and slow 
motion. Clearly this situation is also realized in the mod- 
els of stick-slip motion where (especially in the context of 
earthquakes) there is a strong separation of time scales 
between fast slipping events and slow accumulation of 
stress. It is therefore advantageous to try to apply these 
techniques to the problem of stick-slip motion. 

In this paper we present a study of the Burridge- 
Knopoff slider-block model jll| with the Coulomb fric- 
tion law. We will show that for sufficiently slowly spa- 
tially varying displacement variable the dynamics of the 
system is dominated by self-sustained traveling shock 
waves. We will study the properties of these waves and 
reformulate the dynamics of the system in terms of their 
motion. 

Our paper is organized as follows. In Sec. II we intro- 
duce the governing equations for the model we study and 
discuss the features of the friction law used, in Sec. Ill 
we construct the solutions in the form of self-sustained 
traveling shock waves, in Sec. IV we reformulate the 
dynamics of the system in terms of the motion of these 
shock waves and study general properties of the reduced 
problem, in Sec. V we analyze two different types of 
solutions, and in Sec. VI we draw conclusions. 



1 



II. MODEL 

The Burridge-Knopoff model consists of a one- 
dimensional array of blocks of mass m resting on a fric- 
tional surface . The blocks are connected together by 
springs with spring constant k c and pulled by a loader 
plate moving with constant speed V via another set of 
springs having spring constants k p (see Fig. [j]). Let us 
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FIG. 1. The Burridge-Knopoff model 

measure the displacement Xi of the «-th block relative to 
the point of attachment of the i-th loader spring. In this 
case the total force Fi acting on the i-th block is given 
by 



Fi — k c (Xi+i + Xi-i — 2Xi 



kpXi 



fii 



(1) 



where fi is the force of friction. The dynamics of the sys- 
tem is completely determined by the equation of motion 
mXi — Fi, provided that the friction law is specified. 
Note that the friction force /, is the only nonlinearity 
in the equation of motion. Clearly the dynamics of the 
system will significantly depend on the particular choice 
of the friction law. Recently, a lot of results were pre- 
sented on the dynamics of the Burridge-Knopoff model in 
the case of velocity- weakening friction law P, pT| , ^2] -^6| . 
A characteristic feature of the Burridge-Knopoff model 
with this form of friction is its highly chaotic dynam- 
ics that occurs on all length scales down to the smallest 
length scale a and therefore, the absence of the proper 
continuum limit [ fTlfl . 

In contrast to most previous studies, here we adopt 
the Coulomb friction law, which is applicable to clean 
dry surfaces (see, for example, p3j). Namely, we will 
characterize the friction force by the value of the static 
friction f r and the kinematic friction f s < f r , where f r 
and f s are positive constants. Thus, the block will remain 
at rest if Fi < f r . When Fi reaches f r the block starts to 
move (slips) and the friction force drops instantaneously 
from the value of f r to the value of f s . When the block 
comes to rest (sticks) , static friction turns on again. Also, 
we will supplement the friction force with the viscous 
friction term /viscous(«) = ctv, where v is the velocity of 
the block and a is a constant (Fig. ||). Note that this 
kind of friction law was considered already in the original 
paper of Burridge and Knopoff and is observed in the 
friction experiments and their analogs [f2l|-|25| . Also note 
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FIG. 2. Friction force as a function of the velocity 

that this model was recently studied numerically in Ref. 
|l7f (although for the parameters different from those 
considered in this paper). 

We would like to emphasize that in contrast to velocity- 
weakening friction law, for which the kinematic friction 
force is a continuous function of the block velocity at 
v = 0, in our friction law the force / has a discontinu- 
ity at v = 0, reflecting the fact that the coefficient of 
static friction is typically larger than that of kinematic 
friction. As a result of this difference, the accumulated 
stress accelerates the block at the onset of motion, releas- 
ing potential energy. As we will see below, this provides 
the sustaining force for the traveling waves studied in 
the following section. Also, as we will show below, this 
ensures the existence of the proper continuum limit as 
a^O. 

Let us now formulate the equations of motion for the 
blocks in the case of the friction law introduced above. 
Recall that the displacements Xi are measured relative to 
the loader plate, so they are defined in the reference frame 
moving relative to the surface with velocity V. Therefore, 
when the blocks are at rest, their equation of motion 
becomes 

Xi = -V. (2) 
When the block slides, its equation of motion changes to 
ml, = k c (X l+1 + AV! - 2X t ) - k p X t - f s - a(X t + V). 

(3) 

The transition to sliding (slip) occurs when the force F% 
reaches the value of f r : 



k c (Xi + i + Xi-i — 2Xi) — kpXi — f r . 



(4) 



In other words, the equation of motion of a block changes 
from Eq. (||) to Eq. (]|) when the condition in Eq. (H) is 
satisfied. Of course, the block comes back to rest when 
Xi = — V during sliding. 

Let us introduce the following dimcnsionlcss quantities 



a V k r V m 



kpXi + f s 

u * = — /WT' (5) 



2 



where a is the distance between the attachment points of 
the loader springs (Fig. [j]). The dimensionless displace- 
ment variable it, is chosen in such a way that the system 
is excitable only at m > [see Eq. (|l|), in order for slip 
to occur, the force Fi must exceed the value f s of sliding 
friction] , while the blocks will always slide for it, > 1 [see 
Eq. (|)]. 

If k c k p , one can naturally go to the continuum limit 
by making a substitution X i+ i + — 2JQ — > a 2 ^^-. 
This is a good approximation when the length scale of 
variation of X» is of order 1 in the new variables. Using 
the variables in Eq. (||), we rewrite Eqs. (^) and (||), 
respectively, as follows: 



u t = v, 

u tt = u xx -u - 2j(u t - v), 
where we introduced dimensionless constants 



Q 



2^/kpT 



Vy/k~n 

fr ~ /s 



(6) 
(7) 



(8) 



and dropped the primes for simplicity of notation. Note 
that in the continuum limit in addition to tracking the 
motion of individual blocks, one also has to follow the 
motion of the slip points. Similarly, one should keep track 
of the positions of the points at which the blocks come to 
rest. The latter are determined by the condition in Eq. 
(|J) during sliding. 

Determining the position of the slip points in the con- 
tinuum limit turns out to be a rather complicated prob- 
lem, since the position of the slip will still depend on the 
local dynamics of the blocks. In the discretized form, 
the slip condition [Eq. (Q)] in the variables of Eq. (^) 
becomes 



H+l 



2u,, 



(Ax) 2 



1, 



(9) 



where Ax = \Jk p jk c . We will get back to this problem 
in the following section where we will show that in the 
continuum limit the dynamics of the slip becomes inde- 
pendent of Ax. 

Equations @) - (|J) are the constitutive equations that 
will be studied in the rest of the paper. As can be seen 
from Eq. (§), 

we chose time and length scales so that 
the speed of sound in the system is equal to 1. The time 
scale is determined by the period of oscillations of an iso- 
lated block due to the loader spring. Note that in the 
continuum limit the characteristic speed of a block dur- 
ing sliding is much smaller than the speed of sound. This 
can be seen from Eq. (|5|) if one assumes itj ~ 1 , measures 
Xi in the units of Eq. (|5|), and uses a natural assump- 
tion that f r — f s <C ak p . The behavior of the system is 
determined only by two parameters: the dimensionless 
dissipation parameter 7 which measures the effect of vis- 
cous friction, and the dimensionless rate of accumulation 
of stress v. In the following we will consider v to be small, 



what expresses the fact that the time scale of accumula- 
tion of stress is much longer than that of the motion of 
an individual block during sliding. 

Finally, let us discuss the applicability of the Burridge- 
Knopoff model in the context of real physical systems 
exhibiting stick-slip motion. In a real system one should 
replace a one-dimensional array of masses between the 
loader plate and the frictional surface by an elastic 
medium of certain thickness h. A straightforward exten- 
sion of the Burridge-Knopoff model would therefore be 
a two-dimensional array of masses connected by springs 
whose one edge is rigidly attached to the loader plate 
and the other slides on the frictional surface. Naturally, 
the thickness of such a medium should greatly exceed the 
microscopic length scale a. The stick-slip motion in such 
a medium is due to surface waves p6| . The dispersion 
of these waves is given by ui 2 — k 2 + (7r/2/i) 2 (l + 2n) 2 , 
where n is an integer and the transverse speed of sound 
was taken to be 1. The main difficulty here is that in- 
stead of a single displacement variable u one has to deal 
with a large number of modes with different n. A simpli- 
fication that is presented by the Burridge-Knopoff model 
consists of lumping up these modes into a single mode 
with an effective stiffness k p [Eq. (0)]. It is clear that the 
dominant modes will be those with small n, so we must 
have k p ~ (a/h) 2 k c <C k c . Thus, in a physically rele- 
vant situation one should consider the continuum limit 
Ax — > of the Burridge-Knopoff model. In short, this is 
merely a simplification of the mathematical handling of 
the elastic dynamics. A much more important assump- 
tion is that the friction response to the motion of the 
medium is instantaneous and has zero correlation length. 
It would be incorrect to think of the "size" of the block 
a as an atomic distance. Rather, the value of the pa- 
rameter a should have a magnitude of the correlation 
(memory) length of the cooperative effects responsible 
for the surface friction (see, for example, [^lj and refer- 
ences therein). Then the response of the system can be 
considered to be instantaneous if the characteristic ve- 
locity of the blocks is greater than a/r s ti c k, where r s ti c k 
is the characteristic sticking time. 



III. TRAVELING WAVES 

Let us now demonstrate that Eqs. (||) - (^) admit solu- 
tions in the form of traveling shock waves with constant 
speed c in the limit of vanishingly small v. But before 
we do that, let us see what kind of dynamics one would 
observe if there are no slip events and the initial distri- 
bution of u is taken to be a sufficiently slowly varying 
distribution uq{x) (the latter ensures that no slip events 
will occur at short times). In this case the dynamics is 
trivial: we will have u(x,t) = Uo(x) + vt, as long as no 
slip events occur [see Eq. (0)] . From this one can see that 
the characteristic time scale of accumulation of stress is 
v 3> 1, when v is small. In other words, for w < 1 on 
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the time scale of order 1 one will not see any motion with 
the distribution of u fixed to Uq [x) + C, where C is a con- 
stant. In particular, if one starts with the uniform initial 
conditions, one will have u = u+, where u+ is some con- 
stant less than 1. The other possible situation when the 
dynamics of the system becomes trivial is steady creep, 
when all the blocks move together with the loader plate, 
so we simply have u = 2-fv ~ for u«l [see Eq. (0)]. 

Consider the system with u — u + and v = 0. As was 
noted in the previous section, when u + < slip events 
cannot occur, so if one slightly moves a particular block, 
it will quickly move to a new equilibrium position and 
no significant changes will occur. A different situation is 
realized for < m+ < 1 when the system becomes ex- 
citable. Then, if a single block is slightly moved, it will 
accelerate, releasing the potential energy accumulated in 
u + . As a result of this acceleration, the force acting 
on the adjacent blocks will increase resulting in slipping 
of the adjacent block that was at rest. This avalanche- 
like process will go on, leading to the formation of the 
shock wave that releases the accumulated stress, so that 
u changes to a new lower value of u_. The situation 
here resembles a great deal combustion of safety fuse dis- 
cussed in the Introduction. Thus, for < u + < 1 a small 
perturbation of u may lead to a significant change of the 
state of the system, switching it from u = u + to u = u_. 
As in other excitable systems [jl]-|| , it is natural to expect 
that such a perturbation will transform at long times to a 
shock wave traveling with constant velocity c. Numerical 
simulations of Eqs. (||) - (^) with v — show that this 
indeed happens in our model. 

Since all the "nonlinearity" in the problem is contained 
in the slip condition given by Eq. (^|) , our system is not 
unlike piecewise-linear model used to study propagation 
of nerve impulses p7| . For this reason the solution in 
the form of the traveling wave can be found exactly for 
V = in the continuum limit. For definiteness, we will 
consider the waves traveling from left to right, so that c > 
0. Because of the reflection symmetry, for each solution 
traveling with speed c there is also a solution traveling 
with speed — c. 

It is clear that the speed c of the traveling wave should 
be determined by the motion of the slip point in front of 
the wave. Let the slip event occur for block s at t = t s . 
Then, for i > s we have Ui = u + , so the slip condition 
from Eq. becomes 



v {t s )=u + -{l-u+){Axf 



(10) 



Since at the onset of the slip u s = u + and du s /dt = 0, at 
time t we will have u s = u+ + 0((t — i s ) 2 ), with u s < u+, 
since right after the slip the block is accelerated by an 
excess force. This means that at some time t s +i the slip 
condition from Eq. (^) will be satisfied for the (s + 1)- 
th block. According to Eq. (|io|), this will happen when 
At = t s+ i — t s = O(Ax). Therefore, on the time scale 
of order 1 this will correspond to the motion of the slip 
point with the speed c = Ax/ At = O(l). 



FIG. 3. The dependence c(u+). Results of the numerical 
solution of Eqs. (A.;) - (|A7| ). Dashed line shows the result of 
approximation by Eq. (JXW- 

To actually calculate the speed c of the front, one needs 
to solve the system of coupled equations of motion for 
the blocks behind the s-th block. This problem in the 
limit Ax — > is considered in Appendix A. The analy- 
sis of this problem shows that the speed c of the wave 
is uniquely determined by the value of u = u + in front 
of the wave. Note that similar situation takes place in 
the models of combustion (see, for example, |18|) and in 
reaction-diffusion systems (see, for example, Pppl). The 
dependence c(u+) found numerically is shown in Fig. [|. 
We would like to emphasize that upon the decrease of 
Ax the speed c becomes independent of Ax, providing a 
well-defined continuum limit for slip propagation. These 
results are also supported by the direct numerical simu- 
lations of Eqs. (J6j> - (0). 

According to Fig. (H), traveling wave solution exists 
for all < u + < 1. The speed c is always greater than 1, 
that is, the considered shock waves are supersonic. The 
latter is also observed in other models of stick-slip motion 
Jl5| , ^6| . Note, however, that the speed c is the speed of 
propagation of the slip and is not related with the actual 
speed of individual blocks in the wave. 

As can be seen from Fig. (||), the speed of the travel- 
ing wave diverges as u+ approaches 1. One can get an 
analytical handle on the dependence c(it+) by expand- 
ing it in the inverse powers of c (see Appendix A). As a 
result, we get the following interpolation formula which 
implicitly determines c as a function of u+ 



u + 



1 



4c 4 



8c 5 Vc 3 ' 



1 



(11) 



This equation gives the dependence c(m+) with a few per 
cent accuracy. 

Let us introduce self-similar variable z = x — ct, where 
c = c(u + ) (Fig. |J). Since the problem possesses transla- 
tional invariance, we can choose the position of the slip 
point to be at z = 0. Similarly, the stick point z = — w, 
with w > 0, will also travel with speed c behind the wave. 
Behind the slip the blocks slide according to Eq. (0), so 
for the wave with speed c we will have 



(c 2 — l)u zz — 2cju z + u = 0. 



(12) 
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The slip condition gives the following initial conditions 
for this equation: 



u(0) = u+, «»(0) = 0. 



(13) 



The analysis of Eq. (|T^) (see Appendix B) shows that the 
structure of the solution changes qualitatively depending 
on whether the value of u + is smaller or larger than the 
critical value u c > 0, where the value of u c is implicitly 
determined by 
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Vc 2 K)-i 

c(u c ) 



(14) 



If we have m+ < u c , the distribution of u will asymp- 
totically approach u_ = at z = — oo behind the wave 
[Fig. |(a)] . 

In other words, the blocks behind the wave never come 
to rest. Therefore, upon passing of the wave the system 
should develop a steady creep for small but finite v. Note 
that this will always happen when 7 > 1 since in any 
case u + < 1, so for these values of 7 the propagation of 
multiple shock waves becomes impossible (see Sec. V). 

A different situation takes place for u + > u c . Then the 
solution has the form shown in Fig. ^(b), and the width 
w of the traveling wave becomes finite (see Appendix B) 



7T(C 2 - 1) 



v /c 2 (l- 7 2 )-l 



(15) 



Thus, in this case the wave propagates as a "self-healing" 
pulse (compare with p8|). Note that for 7 > we have 
c > c(u c ) > 1, so the value of w is bounded from below 
and the width of the wave is always greater or of order 1. 
This justifies the use of the continuum limit for k c 3> k p 
and 7 ~ 1. Furthermore, when u + becomes close to 1, 
the propagation speed c and therefore the width of the 
wave grow, so in this case the continuum limit is justified 
even for k c ~ k p . On the other hand, the width w goes 
to zero when u + becomes small when 7 « 1 [see Eq. 
(|l5|)], so for sufficiently small u + and 7 the continuum 
limit will become invalid for the description of traveling 
waves. 

In the case u+ > u c the value of u = u_ behind the 
wave will be (see Appendix B) 



-u+ exp 



7T7C 



Vc 2 (l-7 2 )-l 



(16) 



From this equation one can see that it_ < 0, so right 
behind the wave the system is in the state in which no 
waves can be further excited. Also, note that when u + 
approaches the value of u Cl the value of u_ rapidly ap- 
proaches zero. 

In the context of earthquakes a shock wave should be 
associated with an individual earthquake. The total dis- 
placement A = u + — u_ which occurs upon passing of a 
wave and which determines the magnitude of the earth- 
quake will therefore be 



A = u 4 



1 + exp 



7T7C 



Vc 2 (l-7 2 )-l 



(17) 



IV. FREE BOUNDARY PROBLEM 

So far we were studying traveling wave solutions in the 
limit v = and u — const in front of the wave. These 
solutions are clearly the fast motions in the system since 
they occur on the time scales of order 1. On the other 
hand, although the accumulation of stress occurs slowly 
for v <C 1, it can lead to significant changes in u at long 
times of order v^ 1 ^> 1, so the latter can be associated 
with the slow motions in the system. Therefore, since 
the fast motions result in the large-scale changes in u, 
the slip-stick events can be considered as singular per- 
turbations to the slow motions. 

Similar situation is realized in other excitable media 
where the role of singular perturbation is played by the 
sharp fronts . A powerful tool in describing the dy- 
namics of such systems is reduction to free boundary 
problem (see, for example, [p^| ). In this approach one 
obtains the relationship between the speed of the sharp 
front and slow variables. Using this idea, let us formulate 
the dynamics of our system as a free boundary problem. 

Let us assume that the distribution of u varies on the 
length scale much greater than 1. Then, introducing an 
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adiabatic approximation, we can assume that the speed 
of the traveling wave is given by the dependence c(u+) 
(Fig. ||) in which now, instead of u + one should use an 
instantaneous value of u right in front of the wave. The 
latter will play the role of the slow variable. Thus, the 
motion of the shocks can be described by the variables 
Xi (t) which give the positions of the «-th shock, and the 
variables s,- = ±1 which give the direction of their mo- 
tion. 

On large length scales the shock will contribute a dis- 
continuity to the distribution of u in the limit v — > 0. It 
can be included in Eq. (^3j) describing slow motions in the 
form of a <5-function. Therefore, in the presence of the 
shocks Eq. (ph can be rewritten as 



u t = v - 



CiA(u)6(x - Xi + Osi), 



(18) 



where A is the amplitude of the shock from Eq. (^) eval- 
uated at u+ = u(xi(t),t), Ci is the absolute value of the 
speed of the i-th shock, and the last term in 5-function 
represents the fact that the value of u is evaluated right 
in front of the wave. Equation (|l8|) simply says that at 
the moment t the value of u jumps from u(xi(t),t) at 
x = Xi(t) to the new value U-(u(xi(t),t)) given by Eq. 

©• . 

Having now defined the evolution of the slow variable 
U, we can write down the equation of motion of the shocks 
in terms of it 

Xi SiC-i^ (1^0 

where Ci is the function of u + = u(xi(t),t) (Fig. |J). 

Equations (|l8|) and (|l9|) are the basic equations of the 
free boundary problem describing the dynamics of the 
system for t C 1 and sufficiently slowly varying initial 
conditions, assuming that no creep occurs in the system 
during its evolution. These equations have to be supple- 
mented with the way of dealing with creation and anni- 
hilation of shocks. From the physical considerations it 
is clear that a pair of shocks moving toward each other 
will annihilate. Similarly, when the value of u reaches 1 
[recall that our distributions of u vary sufficiently slowly, 
so one can ignore the variation of u in Eq. (0)], a slip 
event is initiated at the point where u = 1, creating a pair 
of counter-propagating shocks. These features should be 
incorporated into the free boundary problem and will be 
discussed in the following section. Notice that according 
to Eq. (||) the distribution of u for blocks at rest must be 
a continuously differentiablc function in the continuum 
limit, so the speeds Cj of the shocks given by Eq. (|l9| ) 
will also be continuously differentiable functions of time. 

In writing Eq. (18) we assumed that we always have 
u(xi{t),t) > u c , so that the shock waves are always the 
waves of switching between the blocks with u = u + at 
rest in front of the wave to u — u_ at rest behind the 
wave. In other words, we do not consider the possibility 
of the onset of creep behind the wave. The analysis of 
creep motion is beyond the scope of the present paper. 
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FIG. 5. Periodic wave train. 

The condition of the absence of creep should in fact be 
satisfied for a wide class of initial conditions. In particu- 
lar, it is easy to see that a wave with it+ > u c will never 
be able to reach the points where u < u c if we have 



< 



c{u c ) ' 



(20) 



provided that the regions with u < u c are initially at 
rest. The condition in Eq. ( pp[ ) simply means that the 
wave with speed c will never catch up with the point at 
which u = u c when u grows according to Eq. (^|). 



V. SPATIOTEMPORAL PATTERNS 

Let us now apply the procedure developed in the pre- 
ceding section to spatiotemporal patterns forming in the 
system under consideration. As in any excitable system 
the basic pattern of this kind is a periodic array 
of shocks (wave train) traveling with constant velocity. 
From the results of the preceding section one can see 
that a periodic wave train should consist of sharp fronts 
in which the value of u jumps from w+ to U- followed by 
refractory regions where u slowly recovers back to u+ , so 
the resulting pattern is saw-tooth (Fig. |s[> - In the refrac- 
tory region u obeys Eq. (^]), which in the frame moving 
together with the wave with speed c becomes u z = —v/c, 
where z = x — ct is the self-similar variable. Therefore, 
the solution for u in the refractory region which properly 
matches with the back of one shock at z — and the 
front of another at z — — L, where L is the period of the 
wave train, has the form u = L^ 1 (z + L)u- — L zu+, 
where 



cA 

V 



(21) 



c = c(u + ) (Fig. ||), and A is given by Eq. (17). Equation 
( pl| ) should in fact be obvious from the purely geometric 
considerations. Thus, there exists a family of nonlinear 
periodic traveling wave solutions whose speed and period 
depend strongly on amplitude. Similarly, the period T of 
this solution in time depends on the amplitude simply as 
T = A/v, so the frequency of the shocks at a particular 
point as the wave train passes is lu ~ A^ 1 . The latter in 



G 




FIG. 6. A source of counter-propagating waves. 



fact represents a simple form of the Gutenberg-Richtcr 
scaling law with scaling exponent 6=1 p9[ ]. In other 
words, this exponent would be observed in the system 
under consideration if its dynamics were dominated by 
periodic wave trains. It is interesting to note that the 
value of b actually observed from the earthquakes gener- 
ally lies in the range 0.8 < b < 1.2 |§. 

The above arguments are justified as long as u+ > 
u c , so that no creep develops behind the traveling front. 
Therefore, according to Eqs. (|l7|), and ([2l|), traveling 
wave trains exist only when L > L m i n , where 



L m in = V 1 c(u c )u c 



(22) 



Let us now consider another kind of spatiotemporal 
patterns that is typical of excitable systems — the source 
of counter- propagating waves |l]-|(| . The existence of such 
a solution is associated with the fact that if there is an in- 
homogeneous distribution of u at rest such that the maxi- 
mum of u is located at x = xo, then at some point in time 
the value of u(xq) may reach 1, so that the blocks will 
become unstable creating a pair of counter-propagating 
shocks (Fig. ||). After these shocks moved away from iro, 
the distribution of u can once again have a maximum at 
x = xo, so after some time the value of u will increase 
until it reaches 1 and the cycle repeats. 

These solutions are indeed realized in our system. 
Close to x — xq we can approximate the distribution 
of u as 



u(x,t) ~ 1 — a(x — xq) 2 



(23) 



[see Eq. ([!])], where we also assumed that the slip occurs 
at t = 0. This time-dependent distribution of u gives 
the values of u in front of the shocks at t > 0, after the 
slip occurred, thus determining the velocity of the shocks. 
The analysis of the free boundary problem in this case 
(Appendix C) then shows that right after the slip the 
positions of the shocks will be given by 



xq ± Vbt, 



V + y/v 2 + i 
2a 



(24) 



where the plus and the minus correspond to the shocks 
propagating to the right and to the left, respectively. No- 




a/v 2 

FIG. 7. Dependence a' (a) from Eq. at 7 = 0.3 



tice that Eq. (|24|) in general shows how to treat the sin- 
gularity in the free boundary problem (Sec. IV) associ- 
ated with the creation of a a pair of counter-propagating 
shocks. 

The solution for Xi(t) in Eq. ( pi]) then allows to cal- 
culate the distribution v! after the shocks have passed 



u'(x, t) ~ W-(l) — a!{x — xo) 2 + vt, 
where a' is given by (see Appendix C) 

. 2av(l + k) 

a = —an-\ , 

v + v v 2 + 8a 

u_(l) = - cxp(-7T7/ v /l - 7 2 ) [see Eq. (|l6|)], and 



(25) 



(26) 



7T7 



(1 - 7 2 ) 3 / 2 



exp 



7T7 



r 



(27) 



From Eq. (|27| ) one can see that the value of n lies in the 
range < K < 1. The plot of a' as a function of a at a 
particular value of 7 is presented in Fig. 0. 

From Eq. (Eq) one can see that there is a maximum 
of v! at x = (a' > 0) when the shocks have passed, if 



a < a„ 



where 



v 2 (1 + k) 
2^2 ' 



(28) 



is 



satisfied, after time T = 
the distribution of u will 



If this condition 

w _1 + exp(-7T7/ y/l— 7 2 )) 

look exactly like the one in Eq. ( [2 3D which we started 
with, but with a different value of a. Therefore, Eq. ( |2^ ) 
defines an iterative map for the value of a correspond- 
ing to the solutions for u at times nT, where n is the 
number of iterations. Physically, a source with period T 
will form at x — creating traveling waves propagating 
away from it. Note, however, that even though the pe- 
riod of the source is a constant that does not depend on 
the initial conditions, the solution represents an aperi- 
odic process since after each cycle the value of a changes. 
The latter is represented in Fig. A Observe that we al- 
ways have a' < a. From Fig. [7 one can also see that 
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after many periods the value of a will tend to zero. The 
analysis of Eq. ((2^) (see also Fig. shows that zero is 
the only fixed point of the map a — > a' for any 7. For 
small enough values of a Eq. (|26|) can be approximated as 
a' ~ a[l — 2a(l + K,)/v 2 ]. It is then easy to show that after 
n ^> 1 iterations we will have a n ~ v 2 /[2n(l + k)], where 
a n is the result of n iterations of the map. Note that 
because of the fact that a becomes smaller and smaller 
with each cycle, the distribution of u becomes flatter and 
flatter [see Eq. (|23|)], what means that after a long time 
the regions in the neighborhood of the source will tend 
to synchronize. The characteristic size R of a region in 
which the synchronization takes place can be estimated 
as R ~ an 1 ^ 2 , so we will have R ~ t 1 / 2 . Also note that 
for the same reason if there are two sources with unequal 
values of a < a max some distance R apart, after time 
t ~ R 2 the one with the smaller value of a will overwhelm 
the one with the larger value of a. It is then natural to 
expect that if the initial distribution of u has all maxima 
with a < a max , the system will eventually synchronize 
into uniform oscillations in which all the blocks slide at 



once. 



Different situation is realized when a > a„ 



In this 



case the distribution of u forming after the shocks moved 
away from the origin will be a minimum rather than a 
maximum of u, so the next creation of a pair of shocks 
will not occur at x — xq, but rather at different points 
in space. Therefore, it is natural to expect that for suf- 
ficiently random initial conditions one would see shocks 
created and annihilated at random points in space, cre- 
ating a kind of spatiotemporal chaos. 



VI. CONCLUSION 

Let us now summarize the results of our analysis and 
discuss the relationship of our results with those obtained 
for other models of stick-slip motion. In the present 
paper we analytically investigated the Burridge-Knopoff 
slider-block model with Coulomb friction law, which is 
different from the one used in most previous studies 
P, ^0|p^ |f6[] . This difference is expressed in the fact that 
our friction is a discontinuous function when the velocity 
of the block becomes zero. As a result, the system is ca- 
pable of propagating ultrasonic traveling solitary waves 
in the limit of zero loader plate velocity for any value 
of accumulated stress above the excitability threshold 
u = 0. The existence of such solutions is a character- 
istic feature of excitable systems . 

In the continuum limit the solution in the form of the 
traveling shock wave is independent of the short-scale be- 
havior of the system, except for the precise form of the 
dependence of the wave's velocity on the amount of the 
accumulated stress in front of the wave. The latter in 
fact only weakly depends on the dynamics of the model 
at short length scales if the value of u is sufficiently close 
to the slipping threshold u = 1. Note that one might 



want to calculate the dependence of the wave velocity 
on the amount of the accumulated stress ii+ by formally 
combining Eq. ( p"2| ) at z = with the continuum version 
of Eq. at z — 0. This would give an incorrect result 
c = — u+. This means that although the contin- 

uum approximation may be valid for finding the wave 
profile behind the slip point, one still needs to solve the 
discrete problem at the tip of the wave in order to find 
its velocity. 

In the limit of vanishing loader plate velocity the only 
parameter in the model is the dimcnsionless coefficient 
of viscous friction 7. All our analysis was performed for 
arbitrary values of 7, so it is also applicable to the case 
7 = 0, i.e., when the viscous friction is absent. This case, 
however, has one special feature. The thing is that for 
7 = we have u c — [see Eq. jf^)], so as u approaches 
zero the speed of the wave approaches 1 and the width 
w goes to zero [see Eq. (Eq)]. Therefore, for a fixed 
value of Ax <C 1 the continuum limit will no longer be 
justified for the description of the waves when u becomes 
sufficiently close to zero. Also, when 7 = the amplitude 
A = 2u + of these waves [see Eq. (p"7|)] will be able to 
become arbitrarily small. 

According to the results of our analysis, the solution 
in the form of the traveling shock wave, as well as a pe- 
riodic wave train of a given period, is unique in the limit 
of vanishingly small loader plate velocity, and its speed 
is uniquely determined by the value of u in the front of 
the wave. This is in contrast with the results of Langer 
and Tang and Myers and Langer [|1| for the models 
with velocity-weakening friction where they find a con- 
tinuous family of such solutions and therefore a selection 
problem. 

Note that the systems with velocity- weakening friction 
transform to the model considered by us in the limit of 
small Vf , where Vf is the characteristic velocity of vari- 
ation of the friction force. Indeed, if the characteristic 
speed of a block during the slip exceeds Vf, the friction 
force will drop very rapidly, acting in the same way as 
in our model. To obtain the criterion of smallness of 
Vf, we note that the effect of the velocity dependence of 
the friction is most pronounced in the slip region, where, 
according to Eq. ( |l0| ) , u t ~ Ax — \Jk v jk Cl where we as- 
sumed c ~ 1 , so u is not near the threshold value u = I. 
In the original variables dXi/dt ~ (/ r — f s )/^/k c m, so 
the condition Vf <C dXi/dt gives 



fr fs 

\Jk c m 



(29) 



This formula should in fact be obvious from the physical 
considerations. Similarly, one should recover our results 
for the creep-slip models considered in Q if the charac- 
teristic velocity of variation of the friction force is small 
enough. 

In contrast, if the condition opposite to the one in 
Eq. ( p9| ) holds, no traveling shock waves can be real- 
ized in the continuum limit. This can be seen from the 
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following argument. In a traveling wave the profile of u 
should be described by Eq. ([l2|) in which, in the case of 
the velocity-weakening friction one should drop the term 
—2cju z and replace u by u — 1 close to the slip point. 
There we must have u zz < 0, so, according to the mod- 
ified Eq. (|l2|), we will have c < 1 in the traveling wave. 
On the other hand, the propagation of subsonic waves in 
the system is impossible since from the dispersion rela- 
tion lu 2 = 1 + k 2 for the waves in the absence of friction 
we get that their phase velocity c = ut/k = \/l + k~ 2 > 1 
for all k. On the other hand, numerical simulations show 
that in the discrete models with velocity-weakening fric- 
tion traveling shock waves arc indeed observed. However, 
for Aa; <C 1 these waves move with speeds very close to 
the speed of sound and their width is just a few lattice 
spaci ngs (u nless one is close to the threshold u = 1, see 
also [[l^,[16|). One should therefore expect that in the con- 
tinuum limit only discontinuous shocks traveling with the 
speed of sound are feasible, in contrast to the situation 
studied in the present paper. This is a crucial distinction 
of the Burridge-Knopoff model with velocity-weakening 
friction and the one with Coulomb friction law (Fig. |^) 
studied by us. 

As was already discussed above, the propagation speed 
of the wave depends on the short-scale dynamics of in- 
dividual blocks. Therefore, in order to be able to per- 
form numerical simulations of the models under consid- 
eration in the continuum limit, one has to use very small 
discretization of time. This should make direct simula- 
tions of our model rather difficult. On the other hand, 
we showed that if the initial conditions vary sufficiently 
slowly, the dynamics of the system reduces to the motion 
of individual shock waves (see Sec. IV). Therefore, by re- 
formulating the dynamics of the problem in terms of the 
free boundary problem one can dramatically increase the 
speed of the simulations, thus being able to observe the 
dynamics for much longer times. This should be impor- 
tant in the statistical analysis of our s yst em , especially 
in the context of earthquakes (see also 

Let us see how the results obtained from the studies 
of our model should translate to real stick-slip problems. 
Consider, for example, a Bristol board of thickness h — 
2 mm lying on the surface made of the same material 
under external pressure p and pulled at the top plj ] . The 
coefficients fi r and fi s of the static and the sliding friction 
for these surfaces were found to be /i r = 0.37 and fi s = 
0.31, respectively \21\. The memory length which we 
should use for the distance a between the blocks was 
found to be a = 1 ^tm. A single block in the Burridge- 
Knopoff model should be associated with a piece of the 
board of dimensions a x a x h, so the mass m of the 
block should be m = 2.4 x 10~ 12 kg, where we used 
the density p — 1.2 g/cm 3 . From the speed of sound 
c = ay/ k c /m [Eq. (0)] we can calculate the values of 
k c and k p = k c (ira/2h) 2 (see the end of Sec. II). Using 
c = 4000 m/s, we find k c = 3.8 x 10 7 N/m, and k p = 
2.3 x 10 1 N/m. The friction force jump per single block is 



/,. — f s = pa 2 {ii r — fj, 8 ). One can see that the condition of 
Eq. (||) with V f = 10" 5 m/s (see fH) is satisfied if p = 
2 x 10 6 Pa = 20 atm, what gives f r - f s = 1.2 x 10~ 7 N. 
For this value of p the characteristic sliding velocity will 
be A ~ (f r — f s )/y/k p m — 1.6 cm/s. The displacement 
AA in a single wave can also be estimated as AA = 
2(/ r — fs)/k p = 1 x 10~ 8 m. The latter quantity turns 
out to be quite small. Note that both AA and A increase 
as p increases. Also, smaller pressures p would be needed 
for smaller values of h. From |21| one can calculate a = 
6.7 x 10 -6 kg/s for this value of p, leading to 7 = 0.45. 

In the present paper we studied perfectly homogeneous 
system without any external noise. It is clear that in real 
systems some degree of noise should be present. It is 
interesting to know what effect the introduction of ran- 
domness will have on the behavior of the traveling waves. 
Here one can distinguish two different situations. If one 
adds some small randomness in the parameters of in- 
dividual blocks, such as the friction coefficients, spring 
constants and so forth, one would expect the renormal- 
ization of the speed of the traveling waves. The other 
kind of noise would result in sudden transitions of the 
blocks at rest to sliding. This kind of noise can be read- 
ily incorporated into our free boundary problem by intro- 
ducing generation of pairs of counter-propagating waves 
at random points in space. Here the open question is 
how the frequency of this generation should depend on 
the distance to the slipping threshold. It is also inter- 
esting how the presence of such a noise would affect the 
wave statistics. 



APPENDIX A: 

Since the variation of Ui leading to the slip event is of 
order Aa; 2 and is small for small Ax, for the blocks near 
the slip point one can neglect both the variation of Ui and 
the term —2jdui/dt in Eq. (Q), so in the limit Aa; — > 
the right-hand side of this equation simply becomes 
u+ — 1 = const. For convenience, let us introduce the 
new variables r and Vi, such that 



t = t s +rAx, Ul = u + - u + (Ax) 2 K + yJ • ( A1 ) 

Since at r = the s-th block just slipped, we must have 
[see Eq. ©] 



v s -i(0) 



1 — Uj 



(A2) 



with the initial conditions 

».(0) - 0, v' s (0) - 0, (A3) 

where the prime denotes the derivative with respect to 
r. 

In the limit Aa; — * Eq. (Q) can be written as 







d 2 Vj 

dr 2 ~~ 
d 2 v l 
dr 2 



-i + Vi-i - 2vi, i < s, 
t 2 

i-l 



(A4) 
(A5) 



These equations essentially describe a linear array of unit 
masses connected by springs with spring constant 1, with 
the s-th mass attached to a rigid wall and acted upon by 
a time-dependent force — r 2 /2. 

If the wave moves with speed c, after time r = 1/c, 
which corresponds to At — Ax/c [see Eq. QA1])], it should 
move the distance Ace to the right. This means that at 
time t s+ i = t s + At we must have Uj_i(f s ) = Ui{t s + At), 
u 'i-i(ts) — u 'i(ts + At), or in terms of the new variables 



_ 1 (0)=u i ( C - 1 ) + - 5) 



1 

2d 2 '' 

v' i _ 1 (0)=v' i (c- 1 ) + -. 



(A6) 
(A7) 



Note that these conditions introduce a feedba ck i nto Eqs. 
([A3) and (|Al). Also, the solution of Eqs. should 
be matched with the continuum profile behind the tip, 
which satisfies Eq. (jl^). We should therefore have 



fc 2 



V s -k = 



2(c 2 



1 



fc»l. 



(A8) 



Equations ( |A3| ) - (A7) can be solved numerically by 
an iterative procedure. Then, calculating the value of 
w s (c~ 1 ), using Eq. (|A6|) a s a function of c, one can relate 
it to u+ through Eq. ( JA 2| ) . The results of this numerical 
solution for c > 1 are presented in Fig. (||). These results 
also show that close to u + = we have 
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(A9) 



The problem in Eqs. (A3) can be treated analytically 
by expanding its solution in the powers of c _1 . Indeed, 
i f c is larg e, fo r r < c" 1 <C 1 the right-hand sides in Eqs. 
( A4) an d (]A5 |) can be neglected, so the solution of Eqs. 
(A3) - (|A7|) for fc > in the leading order in cr 1 will be 



,(°) 



2c 2 



k.T 

c 



(A10) 



Physically, this corresponds to the situation in which the 
springs between the blocks are absent. 

To calculate the first order correction vi , we use r 
in Eq. QA9) to obtain 



(0) 



T 

Ad 2 



6c 



T 

24' 



(AH) 



This solution is then used to fix through Eqs. (A6) and 
(A7) the initial condition in Eq. (A4) with i = s — 1 to 



the next order in c 
„.(°) 



Using the solutions vi ^ , v < f} 1 



and 



2 in Eq. (A4), one can then calculate v 



(i) 



.,(!) 



5r 
6^ 



T 

2d 2 ' 



(A12) 



The knowledge of and vi then allows to calculate 
the next order correction vl by substituting them into 



Eq. (|A5|) 



,(2) 



3r 2 



5r J 



16c 4 36c 3 



T 

360 



Substituting the obtained solution for v s into Eq. ( [A2] ) 
and using Eq. (A6) with i = s — 1, we get 



1 

2^2 



3 

8? 



1 — tl4 



16c 6 



(A14) 



Note that this procedure can be continued to an arbi- 
trary order in c _1 , although the calculation then becomes 
rather cumbersome. Let us only quote the result that the 
correction to the left-hand side of Eq. (A14) due to the 
next order terms will be .2719 /c 8 . 

Finally, combining Eq. (Al4)withEq. jA9|) , we arrive 
at the interpolation formula given by Eq. (Illl) . 



APPENDIX B: 

Here we solve Eq. ( |l2| ) with the initial conditions given 
by Eq. ( |l3| ) for arbitrary 7. Let us introduce the new 
variable £ = z/vc 2 — 1 and a constant 7 = cy/^c 2 — 1, 
where c is a function of u + (Fig. |^). Then, Eq. ( |l2| ) can 
be rewritten as 

Utf - 2tu c + u = 0, (Bl) 

with the same initial conditions written in terms of £: 

u(0) = u+, u ( (0)=0. (B2) 

Two situations are possible here. If we have 7 > 1, 
both roots of the characteristic equation of Eq. (Bl) are 
real, so the solution has the form 



ae 



?(7+V7 2 -l) + & e «7-V7 2 -l). 



(B3) 



Using the initial conditions from Eq. ( B2 ) , we obtain for 
the coefficients a and b 



7 




7 



(B4) 



From the solution we see that behind the wave the distri- 
bution of u asymptotically approaches zero. This means 
that we have U- = 0, but the system comes to rest only 
at z = — 00, so the width w of the wave is infinite. The 
form of the solution is shown in Fig. 0(a). 
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In contrast, when 7 < 1, the roots of the characteristic 
equation of Eq. (Bl) become complex conjugate, so the 
solution has the form 



u = ae^ cos^V 1 ~ 7 2 ) + be^ sm{£y/l - 7 2 ). (B5) 



Substituting the initial conditions from Eq. (B2), we 
obtain 



u+, 



u+1 



;2 



(B6) 



In constructing the full solution one should also take into 
account that the block will stick back when it comes to 
rest, what will ha ppen when up = for vanishing v. 
According to Eqs. (|B5| ) and (B6), this will occur at z = 
z s = —w, where 

(B7) 



^ 2 



^7 



Behind the stick point the value of u will remain constant 
equal to u_, with 



the constant b has to be determined. Substituting this 
expression for Xi (€) into Eq. ( |C2[ ) , we obtain the following 
quadratic equation for b: 



ab 2 - vb - 2 = 0. 



(C3) 



According to the definition of b, only the positive solution 
of this equation has physical meaning. Taking this into 
account, we obtain that the value of b is given by Eq. 
(II)- 



From the solution of Eq. ( p2[ ) one can see that the 
shock reaches point x at the moment U(x) = x 2 jb. Right 
after the shocks have passed the new value of u at the 
moment ti{x) is given by Eq. (|l6|). Therefore, according 
to Eq. (|^) the distribution v! at a point x and time 
t > ti(x) behind the shock is given by 



u'(x, t) = u_ ti(x)^j] + vt — vti(x), 



(C4) 



where u is eval uate d right before the shock reached x. 
Since u in Eq. (C4) is close to 1, we can Taylor expand 
the dependence in this equation and keep only 

the first two terms. As a result, we obtain 



-u + exp 



7T) 



7 2 , 



(B8) 



Going back to 7, we can rewrite the equations above as 
Eqs. ( fL5| ) and (fig), respectively. The form of the solution 
in this case is shown in Fig. |J(b). 

The solutions obtained above are the exact traveling 
wave solutions of Eqs. (j6j) — (^|) in the case v = 0. As can 
be seen from the construction, traveling wave solution is 
unique for any given value of u+. 



u' (x,t) ~ w_ (1) + k ( ax 



vx 
~b 



vx 
~b 



vt, (C5) 



where k — —du_/du + evaluated at u + = 1. Using the 
explicit expression for u_ together with Eq. (CI) in Eq. 
(p6|), one arrives at Eq. fl2jj). Then, using the value of b 
from Eq. (p4|) , Eq. (C5) can be transformed to Eq. (I2F 



APPENDIX C: 

When a pair of counter-propagating shocks is created, 
their speed will be large, since in the neighborhood of Xq 
the value of u will be close to 1. This allows us to use 



Eq. (A14) and write 



y/2[l-u( Xi ,t)] 



(CI) 



for Xi close to xq. Note that Eq. ( |Cl[ ) expresses the fact 
that for u close to 1 the coupling of the blocks through the 
longitudinal springs becomes inessential (see appendix 

Without any loss of gene rali ty, we can put x$ = 0. 
According to Eqs. (|l9| ) and (CI), we have approximately 



Xi 



± 



vt) 



(C2) 



for not too large \x%\ and t, with the initial condition 
Xi(0) = 0. The solution Xi(t) that satisfies this equation 
and the initial condition has the form Xi = ±Vw, where 
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